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t^. Abstract: We present a global fit to single- and double- inclusive suppression data of high-p^ 
particles in central Au+Au collisions at top RHIC energy. We also include in this analysis data 
on heavy quarks via their D and B meson semi-leptonic decays (i.e. non-photonic electrons). The 
analysis is based on the parton quenching weights for medium-induced gluon radiation computed in 
the BDMPS approximation then embedded in a hydrodynamical description of the bulk medium. 
Our results indicate that values of the transport coefficient q more than four times larger than 
perturbative estimates are preferred by experimental data. This confirms previous calculations 
based on simpler implementations of the medium geometry or only the single-inclusive suppression. 
We also comment on the statistical compatibility of the heavy quark data within a radiative only 
energy loss scenario, and on the sensitivity of the results to nuclear modification of the parton 
distribution functions (PDFs) and to assumptions on the energy loss during times in the collision 
prior to the hydrodynamical behavior. 
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1. Introduction 

The experimental program at the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National 
Laboratory (BNL) has established the jet quenching phenomenon as a fundamental tool in the study 
of hot matter in heavy-ion collisions |]J - see || for reports on the subject. A variety of observables 
in the large transverse momentum region of the produced particle spectrum is contained under 
this generic name: the suppression of the single-inclusive yields of both light hadrons and non- 
photonic electrons; the suppression of the double-inclusive azimuthal correlations; direct photon- 
hadron correlations; the effects on identified jets, jet shapes; etc. A subset of these observables are 
also available, though with different expected magnitude, at the Super Proton Synchrotron (SPS) 
at CERN ||, and all of these observables will play a major role in the upcoming heavy ion program 
at the Large Hadron Collider (LHC) at CERN. 

Another essential discovery at RHIC is the good description of the data, at small and moderate 
transverse momentum, achieved by a hydrodynamical modeling of the matter produced in heavy-ion 
collisions IJ. In particular, azimuthal asymmetries compatible with a hydrodynamical flow starting 
at very early times and with negligible viscosity have been measured Q. 

Although there is no general consensus, several jet quenching studies || [| |7|, §, |9| require 
an interaction of partons with the medium which exceeds that estimated from the asymptoticaly 
large temperature behavior of perturbative Quantum Chromodynamics (pQCD) |T0| , 11, 12, 13 . 
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Corrections to these estimates (arising from higher orders in perturbation theory) are not known 
for quantities as e.g. q but indications exist that they could be sizable jnj . It is also worth recalling 



here that the relevant scale in the perturbative expansion, both in the medium and in the vacuum, is 
the transverse momentum of the radiated gluon. For an opaque medium case, this quantity acquires 
additive contributions from multiple soft scatterings which can be described in terms of light-like 



Wilson lines 0, [15] . The averages over medium configurations of these non-perturbative quantities 
are usually encoded in the single parameter q. Subleading terms in the high-energy approximation 
involved in calculations of radiative energy loss || are not yet known. Thus, the conclusion of 
the large strength remains, but caution is warranted on the precise interpretation of the actual 
magnitude obtained. 

These findings about the large partonic interaction strength in jet quenching and the hydrody- 
namical behavior of the bulk matter have led to the claim that the produced medium is strongly 
coupled [IB]. Such a conclusion is of critical importance for our understanding of QCD. Therefore 



it must be substantiated by detailed studies which should consider: i) As many experimental ob- 
servables as available; ii) A detailed modeling of the medium compatible with experimental data 
on soft particle production; iii) A statistical analysis of the uncertainties in the constraints on the 
medium coming from both the experimental data and the theoretical implementation of energy loss 
and medium modeling. 

The goal of the present paper is to constrain the properties of the medium which can be 
studied through the description of its effects on high-p-^ particles. The medium is described by a 
hydrodynamical model adjusted to fit experimental data on the soft part of the spectrum. The 
modeling of the energy loss, considered as purely due to medium-induced gluon radiation, is that 



proposed some years ago [17], [18|, g(J and used in several phenomenological studies of RHIC data 
|5|, |5J, [], ^1] . The coupling of a hydrodynamical medium description to an energy loss calculation has 
already several antecedents [22], |23|, [24], [7|, |8|, [|. Fits to experimental data using ^-minimization 
of statistical uncertainties for single- and double-inclusive data were presented in [2B], and including 
statistical and systematic uncertainties in pTL |2T|. 



In this work, we extend previous studies by both including a hydrodynamical modeling of the 
medium and considering the case of inclusive pion, dihadron and heavy quark suppression within 



the same analysis. The energy loss distributions will be taken from EH] for light quarks and gluons 



and from for heavy quarks. The observables related to direct photons-hadron correlations and 



identified jets have not been accessible at present energies until recently [g^, |30|. However, additional 
statistics and refinements in the experimental analysis, and also on the theoretical implementation 

are 



in particular, the development of a Monte-Carlo approach to jet medium-modifications [31 



required before physics conclusions can be drawn. For this reason they will not be considered here. 

The plan of the paper is as follows: In Section |2| the hydrodynamical model of the medium 
will be briefly described (with full details provided in p2| , |33|). In Section |3] the pQCD formalism 
used to compute single- and double-inclusive particle production will be outlined. In Section £| the 
implementation of energy loss will be reviewed. Finally, our results will be shown and discussed in 
Section |5], and the conclusions presented in Section ||. 
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2. Hydrodynamical medium modeling 



The energy loss of an energetic parton is determined by the properties of the medium along the 
parton trajectory ||. In relativistic heavy- ion collisions, the medium itself evolves dynamically 
together with the propagation of high-pr partons. Full three-dimensional (3D) ideal hydrodynamics 
32] , |33| is employed in our calculations to describe the space-time evolution of the hot matter created 



in Au+Au collisions at top RHIC energy (^/snn = 200 GeV). 
We solve the equations of energy-momentum conservation 

dfjT^ = (2.1) 

in full 3D space (r, x, y, rj s ) under the assumption that local thermal equilibrium is reached at 
an initial time To =0.6 fm/c and maintained thereafter until freeze-out. Here r, rj s , x and y are 
proper time, space-time rapidity, and two transverse coordinates perpendicular to the beam axis, 
respectively. Ideal hydrodynamics is characterized by the following energy-momentum tensor: 

= (e + P)uV - Pg^, (2.2) 

where e, P and m m are energy density, pressure and local four velocity, respectively. We neglect 
the finite net-baryon density which is small near mid-rapidity at RHIC where all calculations will 
be performed. For the high temperature (T > T c = 170 MeV) Quark-Gluon Plasma (QGP) phase 
we use the equation of state (EOS) of an ultra-relativistic massless parton gas (u, d, s quarks and 
gluons) with a bag pressure B: 

p = i(e-4B). (2.3) 

The bag constant is tuned to be B^ = 247 MeV to match the pressure of the QGP phase to that 
of a hadron resonance gas at critical temperature T c = 170 MeV. On the other hand, a hadron 
resonance gas model including all hadrons up to the mass of the A(1232) is employed for the low 
temperature phase (T < T c ). The hadron resonance gas EOS employed in this study implements 
chemical freeze-out at T c h = 170 MeV for simultaneous reproduction of particle ratios and 



spectra at low px- We assume that the hadronic matter behaves like a perfect fluid all the way 
down to the kinetic decoupling at T dcc = 100 MeV. 

The initial entropy distribution in the transverse plane is proportional to a linear combination 
of the number density of participants, p part , and that of binary collisions, p co ii |33] as given below: 



s oc (1 - a;)p part + xpcou. (2.4) 

The proportionality constant and the fraction of soft and hard components (specified by x = 0.15) 
are so chosen that the centrality dependence of charged particle multiplicity observed by PHOBOS 
35] is reproduced J34[ . 



For the initial condition of the flow velocity, Bjorken's scaling solution [33], namely u x (tq) = 
u y{ T o) — u vs( T o) — 0, is employed. The hydrodynamic model has been tested against experimental 
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data at RHIC energies such as centrality and dependencies of spectra for pions, kaons and 
protons up to 1.5 GeV/c We will use these hydrodynamic solutions, in particular the 



local energy density e along the trajectory of partons, to fit various high-p^ data and to reliably 
constrain the transport coefficient q. The hydrodynamic results such as temperature and energy 
density as functions of the space-time position are publicly available [[37] . 

3. Perturbative cross sections 

The production of a hadron h at transverse momentum p? and rapidity y can be described by the 
factorized expression 

/ ^-^^^^^^^^//(x^Q 2 )^-^^^^), (3.1) 
J x 2 z z — ' dt 



dpxdy 



where f^(x, Q 2 ) are the parton distribution functions (PDFs) for a colliding object, either a proton 
or a nucleus A, x is the momentum fraction of the hadron carried by the parton, and Dk^h{z, n\) 
are the fragmentation functions for a parton k into a hadron h, with z the momentum fraction 
of the parton carried by the hadron. We utilize CTEQ6L and CTEQ6M free proton parton 



densities depending on whether the calculation is at leading order (LO) or next-to-leading order 
(NLO). We take the factorization and renormalization scales as Q 2 = (pt/z) 2 and the fragmentation 



scale as fip = Pt for light meson production, see e.g. |39[. For the case of heavy quarks, we set all 
scales equal to Q 2 = ji 2 F = M 2 + p\ for consistency with the comparison to experimental data in 



the proton-proton case |40]. When the fragmentation takes place in a medium, the corresponding 



medium-modified fragmentation function is modeled as pi 



&) = J ^P(C)D- h (j^, /4) , (3-2) 

where P(C) is the probability distribution of a fractional energy loss, £ = AE/E, of the fast parton 
in the medium |42| and the vacuum fragmentation function, Dl a f th , is taken from a known set, in 



this case from de Florian, Sassot and Stratmann [33] . Eq. ( |3.2| ) assumes that no modification of the 
virtuality evolution takes place in the medium, the only effect being the energy loss of the leading 
parton when traveling through the medium. This assumption has been studied recently |44|| , mainly 
motivated by the need of a more sophisticated approach for more exclusive observables. However, 
the fact that the medium- induced gluon radiation does not present collinear or soft divergences, 
due to formation time effects, indicates that the corrections to ( pT2| ) are 0(l/px)- 

In the case of proton-proton collisions, the description of the data using Eq. ( |3.1|) is good even 
at lowest order (LO, with appropriate K-factors) for light hadrons produced at central rapidities 



For particles containing heavy quarks, the Fixed-Order plus Next-to-Leading-Log (FONLL) 
approach (|5| is used to resum the potentially large log(Mq/pr) terms, with Mq the mass of the 



heavy quark. This scheme leads to a good description of the Tevatron data [|4J], |47], |48|. At RHIC, 
both the PHENIX and the STAR Collaboration have measured the yield of non-photonic electrons 
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|49| . |50f . which are expected to be dominantly produced by the semileptonic decay of the heavy 
quarks inside the hadrons. In this case, the comparison of data with FONLL is reasonably good for 
PHENIX, while the STAR data appear underestimated by a factor ~ 5. The discrepancy between 
the two experimental data set should obviously be better understood before addressing in more 
detail the agreement with the theoretical prediction. 

The relative contribution of charm and beauty to the final yield of electrons is not known 
precisely, as it suffers from large theoretical uncertainties due to the variation of the actual value 



of the mass of the quarks and unknown higher orders [40], see |5T| for recent experimental efforts. 
This uncertainty is particularly relevant for the medium effects as we will discuss below (also |52|| ). 

Besides single-particle inclusive cross sections, we are interested in the case of two-particle 
correlations. The corresponding expression is similar to Eq. ( |3.1|) , but now two partons created in 
the hard scattering fragment independently. The fragmentation of each parton is again described 
by D™ h or D%^ h for the vacuum and the medium respectively. Notice that here we will be 
interested in two-particle correlations which are well separated in azimuthal angle (back-to-back 
correlations), so that correlations from the fragmentation of a single parton into two final hadrons 
are negligible 



4. Energy loss implementation 

The medium-modified fragmentation functions are computed using the Salgado- Wiedemann quench- 
ing weights P(AE/u c , R) |p0| . They depend on two quantities, u c and R, which for a static medium 



of length L and transport coefficient q, are given by 

u c = ^qL 2 , R = u c L. (4.1) 

In a dynamical medium we make use of the following scaling relations |^, 0, [J 

u e c ff (x ,y , r prod ,0) = J d£ (4.2) 

[gL] e ^(x o ,y o ,Vod,0) = J d£q(£), (4.3) 

Wfft A\ 2[^(x O ,?/o,Tp rod ,0)] 2 

R 11 (x , y , T prod , 0) = 7T7TT777 ~T , ( 4 -4) 

corresponding to an equivalent static scenario to write the effective quantities. These quantities 
are for a parton produced at mid-rapidity where the measurements under consideration have been 
performed, and at transverse position (xo>2/o) and time r prod . Each parton traverses the medium 
with a straight-line trajectory given by azimuthal angle and with local transport coefficient <?(£). 
As such, the problem is reduced to specifying the relation between the local value of the trans- 
port coefficient at a given point, £, on the trajectory and the medium properties as given by our 
hydrodynamical model. Based on the estimates in [10], we define 



q(0=Kq QGP (O^K-2e 3 /\0- (4-5) 
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The energy density along the parton trajectory e(£) is then taken from the hydrodynamical model 
of the medium and K is a fit parameter when comparing to experimental data. This parameter 
then can indicate any departure with respect to the ideal gas <?qgp estimated in flOf . 

The production point of the hard scattering is characterized by a production weight w(x ,y ), 
as computed through the product of the profile functions of the two nuclei, 

w(x ,y ) = T Au (x ,y )T Au (b- (x ,y )), (4.6) 

with the profile function given by a 3-parameter Fermi distribution at a given impact parameter 
b taken from Ref. [56|. The trajectory is then defined by an angle <fi and the average of a given 
quantity O is computed as 

(£>) = — [ d(fidxody w(xo,yo)0(xo,yo,(i)), (4.7) 



with N = 2n j dx dy w(xQ,y ). For example, we may compute the average value of q as follows: 

[[qLYf f{x ,y ,<f>) ] 
\x ,y ,4>) 

The medium-modified fragmentation functions are computed in the same manner: 



(g) = T7 / d(f)dxody w{xo,y ) eff — — . (4.8) 

- /V J [Xn, yn, (D) 



Wfe4)>4 / d<j>dx dy w( X(h yo) f ^P(x ,yoA,()D v k Z h & 



(4.9) 



where P(xo,yo,<p,e) is the corresponding quenching weight for the effective parameters Q4.2p and 
( [4.4[ ). In the case of two-particle correlations, two fragmentation functions are computed for particles 
traveling in opposite directions - i.e. measuring different medium lengths and density profiles: 

DA(x , yo, 0, z, » 2 F ) = J ^P(xo, l/o, 0, QDlZ h (j^, fii\ , (4.10) 
with the average given by 

^ J d4KteQdyQw(xo,y Q )iy^^ (4.12) 
4.1 Energy loss for times prior to hydrodynamic behavior 

In a dynamical medium like the one considered here, there is an ambiguity on the value of the 
transport coefficient, defined by Eq. (4J5), for values smaller than the thermalization time r . One 
extreme case is to take q(£) = for £ < r . The absence of any energy-loss effect for these very 
early times is a strong assumption since thermalization is not necessary in the approach in which 
the quenching weights have been computed @. To quantify this uncertainty, we consider three 
different extrapolations for the time from the hard parton production to the thermalization time 
(r = 0.6 fm/c in the present calculation): 
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• Case i): g(£) = for £ < r ; 

• Case ii): = g(r ) for £ < r ; and 

• Caseiii): q(£) = g(r )/£ 3 / 4 for £ < r . 

These extrapolations range from the most extreme assumption of no effect at all before the thermal- 
ization time (Case i)) to a continuous interaction from the production time (taken to be r pr0( j ~ 0.04 
fm/c) with a free-streaming medium with energy density dropping as e(£) ~ l/£ (Case hi)). As it 
turns out, the extracted value of the parameter K in Eq. ( f4.5p is only modestly sensitive to different 
assumptions about the very early time. 

Both the quenching weights and the hydrodynamical profiles are taken from the publicly avail- 
able codes [^7[ and |37[] , respectively. 



5. Results and discussion 

The experimental data used in our analysis are given in terms of the nuclear modification factors 

dN AA /dp 2 T dy _ D AA (z T ,p^ lg ) 

AA ~ (N coll )dN pp /dp^ AA ~ D pp {z T ,p¥*) [ } 

for single- and double-inclusive measurements, respectively. They are computed using the corre- 
sponding medium-modified fragmentation functions, ( |3.2|) , in the cross sections, ( p.l|) , for the AA 
case and the corresponding unmodified (vacuum) ones for the proton-proton case appearing in the 



denominators of ( |5.1| ). In ( |5.1|) , the hadron-triggered fragmentation function is defined as fl58" 

T h\h,2 trig^jtrig An ,ass< 11 ' A . ,hssoc 

da h A \/ dy^dp T 



fin 1 2 lf1ii trl ^f1r) Z f1ii assoc f1r) a 

D AA (z T ,p^) = tn g ^ /dy dp T dy dp T 

' hx /^,triB^/ trig V ) 



Here. z T = p™ soc /p™ g , and we use as factorization scales the px of the hadrons. We highlight that 
the quantities defined in this manner are, in general, very different from the fragmentation functions 
measured, for example, in e + e~ annihilation at similar virtualities due to the strong bias produced 
by triggering on a steeply falling perturbative particle spectrum. 

The data sets included in our analysis are all for Au+Au collisions at top RHIC energy, ^s^n = 
200 GeV and are listed below. 

1. 'Light': PHENIX data on R AA for vr° with > 5 GeV in the - 5 % centrality class [[ZTJ, 
14 data points. 

2. 'Dihadron': STAR data on I AA for Au+Au/d+Au for charged particles in the 0—5 % centrality 
class ||, for p T rig = 8-15 GeV and p T rig > p^ ssoc > 3.0 GeV, 6 data points 1 . 



1 For computational time reasons, we make the calculations with fixed Py Ig = 10 GeV. We have checked that the 
results are not sensitive to changes in this value. With this choice of Ig and p|? soc , the Daa(zt,Pt 1& ) measured 
f39| in d+Au collisions at top RHIC energy is well described. 
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3. 'Heavy': Raa for non-photonic electrons with > 3 GeV, for — 10 % from PHENIX ||60|| . 
11 data points; and for — 5 % from STAR 9 data points. 

In order to obtain the best value of the i^-parameter (or equivalently q) we perform a modified 
X 2 -analysis that treats statistical and systematic uncertainties in a consistent manner, as detailed 
in [^TJ. For any set of data with central values {?/«}, three types of uncertainties aft = A (statistical 
plus uncorrelated systematics added in quadrature), B (correlated systematics) and C (normaliza- 
tion uncertainties in %), are considered and a modified x 2 for each set of theoretical values {yf 1 } 
(characterized by one given value of parameter K in the present case) is obtained by minimizing 



x 2 



v- [m + e B a f + e c ViO c - yf] 2 



cxp. points 



+ 4 + 4, (5.3) 



* = a? f y^^f+egy^ (5 4) 

with respect to the two parameters Eb, £c- Note that our calculations currently do not explicitly 
include theoretical systematic uncertainties in the above formulation. Instead we consider various 
modifications of the theoretical calculations and their impact of the experimentally constrained 
results. In the case of the the 'heavy' set, there is a theoretical uncertainty of ±8 % coming from 



the variation of the renormalization and factorization scales and the heavy quark mass, see [[52 
Since as will be shown, the theory tends to predict less suppression than the experimental data, we 
take the lowest values inside this band as the theoretical input {yf 1 } for the fits. 



It is notable that in the 'light' case, we follow the treatment of systematic uncertainties in |21 



where the type B uncertainties allow a tilt of the data points (considering cases of correlation or anti- 
correlation as a function of p?)- In the 'dihadron' case, we include the statistical uncertainties quoted 
in and an additional 7% normalization uncertainty as used in the constraint analysis [Z7| . In 
the 'heavy' case, both experiments quote systematic uncertainties which are included in the analysis 
here. Note that we consider the type B uncertainties to move the data points up and down together 
(and do not consider the anti-correlation case here). We discuss these uncertainties for the 'heavy' 
case in more detail later. 

We present the results for different values of the proportionality factor between the transport 
coefficient and the energy density, Eq. ( |4.5| ), K = K'/1A6, where K' = 0.5, 1, 2, 3, ... , 20. For the 
'heavy' case we extended the range up to K'= 40. This grid in K is then used to obtain the best 
X 2 values. 

5.1 Light hadrons 

Our first study corresponds to a global fit of Raa and Iaa within the model and assuming a frozen 
value of q for the very early times — Case ii) previously discussed. Additionally, the EKS98 set |xL 



is taken for the nuclear corrections to the PDFs. This fixes our canonical choices for the calculation. 

In Fig. |l| we plot the curves for different i^-factors for both Raa and Iaa and compare with 
the neutral pion data suppression from PHENIX and the corresponding suppression of the back-to- 
back hadron-triggered fragmentation functions from STAR. Shown in Fig. |2| as solid lines are the 



- 8 - 




Figure 1: Nuclear modification factors Raa for single- inclusive (left) and Iaa for hadron-triggered frag- 
mentation functions (right) for different values of K = K'/1A6, where K' = 0.5,1,2,3,..., 20. The 
theoretical curves are for the canonical calculation (i.e. EKS98 for the nPDF and Case ii) for the treat- 
ment of very early times). Shown with a green line is the curve corresponding to the minimum of the 
common fit to Raa and Iaa data: K = 4.1. 

calculated \ 2 results for both Raa and Iaa comparisons. For our canonical choice, the constraint 
on K agrees within one standard deviation uncertainties. We then sum the two \ 2 values as a 
function of K and obtain a common global fit result of K = 4.1 ± 0.6. 

5.1.1 Comparison of different prescriptions for r < tq 

We now repeat the analysis for the two other assumptions for q for times r < r = 0.6 fm/c. The 
result of this analysis leads to the computed values of the x 2 for light hadrons Raa and Iaa shown 
in Fig. [J (dotted and dotted-dashed lines) and to the corresponding determination of the central 
values and the uncertainty band for A% 2 = 1. The results are relatively insensitive to the choice 
of very early time treatment, except for Raa in Case i). This can be understood by the following 
example contribution: the parton travels for the first 0.6 fm/c without any effect and the strong 
dependence of jet quenching on the in-medium path length forces the value of K (and q) to increase. 

It is notable that in the way the averaged transport coefficient is computed in our approach, 
Eq. (^4.8|) , the actual value of q is different for a given value of K depending on the assumption for 
the initial times in the evolution — this is simply because the average value of the energy density 
is different in the three cases. For example, the values at the minimum K = 4.1 are (q) ~ 6.2, 8.7, 
138.5 GeV 2 /fm for the Cases i), ii) and iii) respectively. This also shows the large ambiguity in the 
determination of an average value of q for an evolving medium: unrealistically large values of (q) 
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Small-T extrapolation 

Case i) 

Case ii) 

Case iii) 



4 6 8 10 

K= q/[2e 3/4 ] 



12 



14 



AA 

Total 



AA 

Total 



Total 



Case i): q(x)=0 for x<X ( 



Case ii): q(x)=q(x () ) for x<x 



Case iii): q(z)=^-j^- for x<x (| 



3 4 5 6 7 8 9 10 11 12 



K= q/[28 3/4 ] 



Figure 2: Left: x 2 -values f° r different values of K for light hadrons and for the three different extrap- 
olations for £ < ro considered here: Case i), Case ii) and Case iii) are plotted with dotted, solid and 
dotted-dashed lines respectively. Red lines correspond to single-inclusive 7r° data from PHENIX (Raa) 
and black ones to the double-inclusive measurements by STAR {Iaa)- Right: the corresponding central 
values (minima of the x 2 ) and the uncertainties computed by considering A% 2 = 1 are shown. 



are possible due to the averaging procedure, in particular when the extrapolation is divergent, as 
in Case iii). This is the reason why in this paper, we have made the choice of quoting only values 
of the i^-factor instead of actual values of the averaged transport coefficients. The fact that our 
results are very solid under the two assumptions ii) and iii) provides strong support to the way in 
which the hydrodynamical medium and the jet quenching calculations are merged. 

5.1.2 The effect of the nuclear modification of parton densities 

We now explore the uncertainty associated with the choice of nuclear parton distribution functions. 
We consider three situations: i) nPDFs as given by the EKS98 analysis [^TJ; ii) nPDFs as given by 



the EPS08 analysis HB2 |; iii) proton PDFs, i.e. without any nuclear modification. These three cases 



comprise two extreme situations 2 . In this study, the impact parameter dependence of the nuclear 
PDF corrections is not taken into account explicitly. This may impact the surface bias from parton 
energy loss, and is in part, why we also consider the no shadowing case as an extreme. 

In Fig. ^ we plot again the value of x 2 for each choice of the nuclear PDFs and the corresponding 
central values and uncertainty bands at the Ax 2 = 1 level. The differences are, in this case, smaller 
than the uncertainties from the modeling of the initial times of medium evolution. Interestingly, 
the results for Iaa are again somewhat less sensitive to these uncertainties than the ones for Raa 
— the reason being the cancellation in the double ratio for the former. 

5.2 Non-photonic electrons and heavy quarks suppression 

We now study the 'heavy' quark case of non-photonic electron data within our canonical set-up, 



2 Notice that for a more formal analysis, the use of the uncertainties computed using the Hessian method as done 
in the recent EPS09 set fl63|] may be used. 
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K= q/[2e V4 ] 




K= q/[2e 3 



Figure 3: Left: x 2 - varues f° r different values of K for light hadrons and for three assumptions on the 



nuclear effects in PDFs: EKS98 (solid lines) |3l]], EPS08 p2[ (dotted lines) and no shadowing corrections 
(dashed-dotted lines). As before, red lines correspond to single-inclusive tt° data from PHENIX (Raa) and 
black ones to the double-inclusive measurements by STAR (Iaa)- Right: the corresponding central values 
(minima of the x 2 ) an d the uncertainties computed by considering A% 2 = 1 are shown. 



i.e. Case ii) for extrapolations at small times and EKS98 for nuclear PDFs. Notice that within 
the radiative energy loss formalism, the energy loss of heavy quarks receives corrections due to the 
mass terms in the propagators whose net effect is a reduction of the amount of radiated energy. It 
is a generic prediction of the formalism that the energy loss of a heavy quark is smaller than that 
of a light quark. The computation of this feature, first pointed-out in Ref. |64|, is performed in 
the approximation M/E <C 1 in Refs. |65, 66|, The quenching weights for the multiple-soft 



scattering approximation used here are found in |57[ . 

It was previously mentioned that the calculation has been performed in the M/E <C 1 limit, 
when M and E are the mass and the energy of the heavy quark. This limit is reasonable for the 
charm quark case in the region under study. The case of the bottom quark is, however, less clear, 
as its mass and energy are comparable for a non-negligible region of the phase space under study. 
Further mass term corrections have not as yet been computed in this formalism. 

We now proceed as in the case of the light quarks and compute, for different values of the 
f^-factor, the Raa of the non-photonic electrons from heavy-meson decays computed within the 
FONLL approach f45f . As expected, the suppression is reduced due to the finite mass of the heavy 
quarks, which is in qualitative contradiction with the experimental data as shown in Fig. As noted 
in Ref. []52 , the unknown relative contribution from charm and bottom decays to the final electrons 
introduces an additional uncertainty in this comparison. Here, an 8% theoretical uncertainty is 
used, as obtained by varying the heavy quark masses and the renormalization and factorization 
scales in the calculation p2| . The theoretical uncertainty band is likely larger, especially for larger 
values of q. This uncertainty is not taken for the case in which only the charm quark contributes 
to the final electron yield as the uncertainty band is much smaller in this case. 

We compute again the x 2 values for different values of the K-f actor. In Fig. we do that for the 
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central theoretical value of Raa diminished by the 8% theoretical uncertainty as stated above. The 
first observation is that the values of x 2 never reach a minimum in the studied region of K (thus 
favoring even larger values beyond the current calculation range). Thus, the 'heavy' data sets appear 
incompatible with the 'light' and 'dihadron' data sets within this theoretical framework. However, 
it is striking that the minimum (within our K =factor range) is Xmm ~ 37 for the PHENIX data 
(with 11 data points) and xLin ~ 1-9 for the STAR data (with 9 data points). 



For the STAR data sample |HJ, around pr ~ 5 GeV, the typical statistical uncertainties are 



of order 32%, the box type systematic uncertainties drawn on the points are of order 29%, and 
the normalization uncertainty is 26%. In the publication, the box type uncertainties are stated as 



uncorrelated [p0[ . If they are point-to-point uncorrelated, we should add them in quadrature with 
the statistical uncertainties (type A). As an exercise, we use the full statistical constraint procedure 
with the theory comparison as a simple second-order polynomial (Raa = a + bxpT + cx pl^). One 
finds a best fit with a = +0.70, b = -0.13 GeV" 1 , and c = +0.01 GeV" 2 with a total x 2 = 0.5. 
This represents an extreme in goodness-of-fit with a p-value = 99.8%. In this paper, we have thus 
treated the STAR uncertainties drawn as boxes around the points as standard type B (assuming 
the points are in fact correlated). Even in this case, the total x 2 = 1-1 remains an extremely good 
match (i.e. reflecting much smaller-than-expected fluctuations of the individual data points). This 
results in a rather improbable situation in that even a value K = 4.1 as determined by the 'light' 
and 'dihadron' data has a reasonable p-value for agreement with the 'heavy' data set discussed here. 
Note that a x 2 test is more stringent in that in order for K = 4.1 to be the true answer, it means 
that K = 30 (for example) is not. Thus, since the x 2 for K = 30 is so small (and much smaller 
than for K = 4.1), it rules out K = 4.1 by the Ax 2 test, but not by a simple p-value check. 

For the PHENIX data sample [|60 |, around px ~ 5 GeV, the typical statistical uncertainties are 
of order 20%, the type B uncertainties of order 19%, and the normalization uncertainty is 7%. In 
this analysis, we consider only the correlated case for the type B uncertainties. We again compare 
the data with a simple second-order polynomial and find the best values as a = +1.6, b = —0.4 
GeV -1 , and c = +0.03 GeV -2 and a total x 2 = 15.2. This is considerably better than the best fit 
with the parton energy loss calculation, but still only has a p-value « 5%. We note that of order 
30% of the total x 2 comes from the single data point at px = 7 GeV. The details discussed above 
regarding the PHENIX and STAR experimental uncertainties may be indicative of not yet fully 
understood systematic uncertainties and/or correlations amongst these uncertainties. Future data 
sets, including with upgraded detectors to specifically tag the displaced vertex from the D and B 
meson decays, will be crucial to resolving these experimental issues. 

In order to further study these data, we make two new computations of the x 2 in Fig. ||] and 
|6]. First we compute the Raa without shadowing corrections. Now the effects are much more 
visible than for light mesons as ant i- shadowing is a 20 — 25% effect here. Second, we compute the 
suppression when only the charm quark contributes to the non-photonic electrons yield, obtaining a 
much better fit as expected — which would be, in fact, compatible with the determination of K from 
light mesons. This second observation provides additional support for the need of an experimental 
identification of the D and B meson suppression. 
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6. Summary and conclusions 



We have performed a fit to several heavy-ion observables in the large transverse momentum region 
of the spectrum. We combine a realistic hydrodynamical medium profile, obtained from fits to the 
soft part of the spectrum, with a jet quenching calculation within the radiative energy loss formalism 
in the BDMPS approximation. The analysis of the light meson data provides the best fit to Raa 
and Iaa with parameter K = 4.1 ± 0.6, for the i^-factor of the fitted transport coefficient with 
respect to a perturbative estimate from Ref. f]10|| . This large value is in qualitative agreement with 



previous estimates, using different medium profiles. This analysis reveals that a common treatment 
of both one- and two-particle inclusive observables leads to statistically compatible results. An 
initial investigation of a sub-set of the theoretical uncertainties reveal the robustness of this fit to 
reasonable changes of these particular assumptions. 

A similar procedure has been applied to study the case of heavy quarks, as measured via non- 
photonic electrons. The baseline spectrum for the heavy quarks has been computed in the FONLL 
approach. Here, the first observation is that the smaller suppression of the formalism, when FONLL 
contributions from charm and bottom are taken into account, leads to a contradiction with the 
qualitative trend of experimental data. The second observation is the large difference of the x 2 
values of the formalism when comparing the theory to the STAR or PHENIX data sets. Values 
of the computed minimum \ 2 are very different with the two experimental data sets, which make 
it challenging to make firm statements about the statistical agreement between theory and data. 
We have also studied the influence of the shadowing corrections in PDFs, and found them to play 
a more relevant role for heavy quarks compared with the case of light mesons. Finally, we have 
studied the case in which all non-photonic electrons come from decays of charmed mesons. In this 
latter case, the comparison gives improved agreement with experimental data due to the smaller 
mass effect of the charm quark with respect to the bottom quark on the radiative energy loss. 

A common fit of all these data ('light', 'dihadron', 'heavy') within a single approach and the 
control over the theoretical and experimental uncertainties remains an open issue in jet quenching 
phenomenology. 
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Figure 4: Nuclear modification factor Raa and the computed values of x 2 for non-photonic electrons from 
heavy quark decays for the same values of the parameter K as in Fig. | — the green line corresponds to the 
same value of K as the green line if Fig. |. The upper panels show the comparison with PHENIX data fl60| | 
and the lower panels the comparison to the STAR data |pC|| . The small difference in centrality, 0-10% in 
PHENIX and 0-5% in STAR is taken into account. The FONLL scheme is used in both calculations. 
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Data: PHENIX 
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Figure 5: Same as Fig. || but without nuclear corrections to the PDFs. The horizontal dashed lines 
in the right panels indicate the minimum of x 2 an d the corresponding values for one standard deviation, 
Ax 2 = 1, and two standard deviations, Ax 2 = 4. Notice that if no local minimum is found in the range of 
K under study (lower panels corresponding to STAR data) these lines are not shown. 
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Data: PHENIX 
charm only FONLL 
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Figure 6: Same as Fig. || but when only the charm quark is taken into account in the calculation. 
Horizontal lines correspond to the minimum of x 2 an d one and two standard deviations as explained in 
Fig. |. 
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